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Abstract 

In this work we detail the application of a fast convolution algorithm computing high 
dimensional integrals to the context of multiplicative noise stochastic processes. The al- 
gorithm provides a numerical solution to the problem of characterizing conditional prob- 
ability density functions at arbitrary time, and we applied it successfully to quadratic 
and piecewise linear diffusion processes. The ability in reproducing statistical features of 
financial return time series, such as thickness of the tails and scaling properties, makes 
this processes appealing for option pricing. Since exact analytical results are missing, 
we exploit the fast convolution as a numerical method alternative to the Monte Carlo 
simulation both in objective and risk neutral settings. In numerical sections we document 
how fast convolution outperforms Monte Carlo both in velocity and efficiency terms. 
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1 Introduction 



Two of the basic problems computational finance has to deal with are the choice of the optimal 
model driving the stochastic evolution of financial variables and, once a good candidate has 
been identified, the search of a reliable way for its fast and accurate simulation. The former 
issue has been widely investigated both by econometricians, mathematicians, and physicists, 



as demonstrated by the increasing literature on this topic, see for example [151 
Tracing back to the work of [2E] and the analysis in [22], empirical studies have shown that 
financial time series exhibit features departing from the Gaussian assumption. In [T7] a de- 
tailed review of the stylized empirical facts emerging in various type of financial markets is 
presented and discussed. These findings are nowadays accepted as universal evidences, shared 
among different markets in different epochs. From the earlier results on cotton prices of Man- 
delbrot or the thick tailed nature of the Dow Jones Industrial Average recognized by Fama, 
very heterogeneous models have been proposed in order to reproduce the degree of asymmetry 
and the excess of kurtosis of the empirical distributions. Approaches directly developing from 
distributional assumptions includes the truncated Levy model discussed in [30J, and those em- 
ploying generalized Student-i and exponential distributions, see [131 [32] . Different mechanisms 
also capturing the observed non trivial structure of higher order correlation functions model 
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the stochastic nature of the return volatility. Continuous time approaches have been exten- 
sively analyzed and range from the fractional Brownian motion, [29], to stochastic volatility 
models, for a review we suggest [21]. Discrete time models include AutoRegressive Conditional 
Heteroskedastic (ARCH) and Generalized ARCH processes, [2DI2], and multifractal ones, [8], 
the latter being inspired by cascades originally introduced by Kolmogorov in the context of 
turbulent flows. Turbulent velocity flows have also led to a series of empirical works testing 
and strongly relying on the Markovian nature of foreign exchange returns, [251 123] . The macro- 
scopic description of the observed phenomena is provided in terms of a Fokker-Planck (FP) 
equation with linear drift and quadratic diffusion coefficients. Processes leading to an equation 
with the same structure characterize several physical systems, as reviewed in [TTJ. Also the 
statistical feedback mechanism proposed in [5J [7] can be recast in terms of non linear diffusion, 
as originally remarked by [33J, who also pointed out potential problems arising when comput- 
ing expectations under power law tailed distributions. Even though these drawbacks have been 
later amended in [38], [34J propose to switch to exponential tailed PDFs. In particular, they 
develop a general approach to generate a Markovian process obeying scaling relations starting 
from a driftless stochastic differential equation (SDE). In the current paper we focus on the 
numerical characterization of these latter processes, of the above mentioned quadratic diffu- 
sion ones, and in their application in financial modeling. To this respect, especially for pricing 
purposes, we need to reconstruct the conditional probability density function (PDF) describing 
the stochastic dynamics in order to evaluate expectations of future payoffs. Yet a closed form 
expression for the density is rarely available. For this reason several numerical procedures have 
been developed and have become the common practice, e.g. binomial and multinomial lattice 
algorithms, Monte Carlo (MC) simulation, and partial differential equation solvers (for a review 
see [14J). We decide to investigate and widely exploit the fast convolution algorithm (FCA) 
introduced in [20] ■ The algorithm applies to Markovian stochastic processes: the repeated ap- 
plication of the Chapman-Kolmogorov equation, and a clever problem re-formulation allow to 
rewrite functional integrals in terms of Fourier and anti-Fourier transforms of the state vector. 
Performing these operations via fast Fourier transform (FFT) algorithm, numerical efficiency 
is achieved and computational complexity is notably reduced. 

The structure of the paper is the following. After introducing stochastic models we have 
chosen to investigate, we detail step by step FCA; in paragraph [2]2 we test its numerical perfor- 



mances against the standard MC approach for different specification of models and parameter 
values. Section[3]is dedicated to financial applications in the context of option pricing. In 3^ 



we 



derive the risk neutral measure for the piecewise diffusion process of |32j , and the exact formula 
for Plain Vanilla pricing. We then consider geometric Asian options by thoroughly develop 



the two dimensional setting required by fast convolution, see paragraph 3.2 We exploit the 



formal analogy between the latter case and the framework discussed by [2S1 El E] to price Plain 



Vanilla options in paragraph 3J3 and in the final part we collect numerical distributions and 
implied volatility surfaces we have obtained to prove the reliability of FCA. We draw relevant 
conclusions and possible perspectives in section [4} 

2 Stochastic processes with multiplicative noise 

Multiplicative stochastic processes we investigate in this work correspond to the class of quadratic 
diffusion described in [TTl [19] and of piecewise linear diffusion introduced in [32] and later on 
rediscussed by pQ in a slightly different flavour. 

The SDE describing the quadratic diffusion dynamics under Ito prescription can be written 



as 
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with X to = initial time condition; dW t is the standard Brownian increment, l/g(t) and e(t) 
are non-negative smooth functions for t > t . We require D 2 = 4 c e(t) — d 2 > with c > 0. 
Equation (fTl) governs the dynamics of a variety of complex phenomena as reviewed in [H] , where 
an exhaustive description of the process in terms of its moments is also provided. For instance, 
when e(t) = e is constant, it is possible to characterize analytically the time-scaling adjustment 
of the process X t , regulated by the choice of g{t), and its convergence to the stationary state. 
If a is non negative or if e is time dependent, X t lacks stationarity; however, in the final section 
of this paper we will discuss an application of this latter case to the context of financial option 
pricing, tracing back to [7j , and later revised by [38] . 

The quadratic diffusion process (fil) can be formally manipulated to reduce it in a more 
convenient form by means of the Lamperti transform, as we will see in Section 2.1 Here, we 
slightly simplify it introducing a new variable beating the time r(t) = f ds/g(s), which we will 
refer to as the integral time. In this new setting the process X T is described by the following 
dynamics 

dX T = (aX T + b)dr + yJcX* + dX T + e(r)dW T , (2) 

with Xq = and e(r) = e(t(r) 

As far as piecewise linear diffusion is concerned, the main property we are interested in is 
scaling, relating returns over different sampling intervals. More precisely, whenever returns are 
rescaled by a factor t H , the shape of their distribution scales according to 
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where G is the so called scaling function. Following [TJ, piecewise diffusion process can be 
defined by means of the driftless SDE 
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Through the analysis of the FP equation satisfied by P(x,t), it can be readily shown that 1/2 
is the only value of H consistent with the scaling assumption (J3J). In addition to this the FP 
equation admits the following solution 
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where a = l/(er 2 e 2 ) and r[a, z] = f^° s a ~ 1 e~ s ds. The interest in previous density is twofold, 
both because of the emergence of a scaling exponent 1/2, in agreement with the empirical 
findings, and because of the evidence of leptokurtosis. [321 EI] derive closed-form option pricing 
formulae under the distribution (p)]); yet they show that consistency between scaling, exponential 
PDF and martingale option pricing requires the replacement of the Ito correction for X t under 
the risk neutral measure with a constant, see section 2 in [31]. We argue this approximation 
is questionable, and we want to perform option pricing not relying on it, therefore we need a 
numerical methodology whose efficiency and flexibility promise to compensate for the absence 
of closed-form solution. 

x By virtue of the properties of <?, r is a monotonously increasing function of t, implying the well-definiteness 
of the inverse function i(r). 



2.1 Fast convolution algorithm 

In this section we review the fast convolution algorithm proposed in |20j . 

Let us consider the generic process X T , whose dynamics is described by the following general 
SDE 

dX T = M X (X T , r)dr + D X (X T , r)dW T , X T=0 = X . 

We start by transforming the process X T into one with unitary diffusion coefficient. This is 
performed via the Lamperti transform, see section 1.11.4 in [26J, defined as 

Z t (X t ,t) ' 
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Under suitable regularity condition Ito Lemma can be applied to Z t (X t ,t) and its dynamics 
turns out to be 

dZ T = M z (Z T ,T)dT + dW T} Z T=0 = 0, (6) 

with 
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where M x is the function M X (X T , r) evaluated in X(Z T ), and analogously for D x . Our aim is to 
provide an approximate expression for the transition probability density function p(z T , t\zq, 0)R 
We introduce an equally spaced time grid = r°, r 1 , . . . , r n = r, with r % = iAr, in a similar 
spirit to the path integral approach, see jTHl E7J H2J |2]. The repeated use of the Chapman- 
Kolmogorov equation in this discrete setting allows to write the transition probability for a 
generic r > as a finite high dimensional integral 

p(z n \z°)~ / ... / Wdz 1 -K{z n \z n ^{z n - 1 \z n - 2 )...'K{z l \z () ), (7) 



n J z n ~ 1 



where z % = z(t 1 ), and it is the short time transition PDF that we chose equal to the Normal 
density 
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By means of the new variables £* = z 1 + Mz(z l , r')Ar, the transition becomes symmetric under 
the exchange of z l+l with £ 4 , i.e. 7t(z l+1 \z l (C, 1 ')) = n ((z l+1 — ^ l ) 2 ). For each one dimensional 
integration appearing in equation Q, we have 

p(z i+1 \z°) = J dz l tt^Iz'M^V) = / df|l^ {{z %+1 - €?) P{z\C)\z°) , (8) 

where p(z l (C, l )\z°) is the density p(z l \z°) evaluated in 2*(£*), and similarly for n. If we introduce 
a numerical integration grid of equally spaced points zj = ^j = z min + jAz for all i — 0, . . . , n 
and j = 0, . . . ,m — 1, where neither z m i n nor A2; depend on the time label i, then the PDF 
^ c ( z P' 1 \Cj) associated to transition of moving from point £j at time r % to point Zji at time r' l+l 



2 From now on we will drop the explicit dependence on the time variable r. 
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is therefore a symmetric Toeplitz matrix 11^ = n|i_j|, with no dependence on the time variable. 
Letting 
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equation ([8| can be approximated as 
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The entries of P are computed by means of the linear interpolation operator 
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becomes 
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(12) 



Matrix multiplications in previous equation are extremely time consuming. Indeed, while mul- 
tiplying a m-vector by the m x m diagonal matrix J requires m operations, and analogously 
for the I operator, multiplication of the II matrix by a vector requires m 2 operations. On top 
of this, the procedure must be repeated at each time step. As a consequence, the dominant 
contribution grows as n x m 2 , and choosing a rather thick grid, computational times rapidly 
explodes. However, the multiplication of a Toeplitz matrix by a vector can be efficiently per- 
formed exploiting algorithms coming from digital signal processing. By embedding II into a 
circulant matrix of dimensions 2m x 2m 
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(13) 



the product of II with a generic vector v is equal to the first m components of Cv e , v c G M 2m , 
v c = (v*, 0, . . . , 0) . Every circulant matrix can be expressed as C = UAU*, where U* denotes 
conjugate transpose of U, whose columns are IP = (l,e _7ry//m , . . . ,e~ m ^ 2m ~ 1 ^ m ) t /\/2m for 
j = 0, . . . , 2m — 1, and A = diag(Co), with Co first row of the circulant matrix. Thanks to this 
result, the product Cv e can be performed exploiting fast Fourier transform (FFT) algorithm 



Cv c = Re [J" 1 (F(Co) • 5F(v e ))] 



where 5F and 3 r_1 are the FFT and anti-FFT operator, respectively, while • is the compo- 
nent wise product. With the adoption of this approach, the computational time is noticeably 
reduced: each FFT computation requires 0(m x log 2 (2m)) operations. Finally, in order to 



compute equation (12) we need to repeat the algorithm at each time step, and on the whole the 
computational burden can be estimated to be of order 0(n x m x log 2 m), which is definitely 
a satisfactory improvement with respect to the non-FFT based procedure. 



2.2 Fast convolution at work: numerical result 

Equipped with FCA we are now ready to approach multiplicative processes previously described. 
We want to check that numerical results obtained by fast convolution converge to the analytical 
solution, when available, or to the PDF reconstructed by means of MC simulation. Moreover, 
we show that FCA provides an estimate of the distribution shape even in those low probability 
regions, such as the tails, which are inefficiently sampled by MC approach. 
Lamperti transform for process p| can be explicitly computed as 
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with A\ = (Ae T c - d 2 )/(4c 2 ), and C° = asinh[(X + d/(2c))/y/A* 
The related drift function is 
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where Xt = (Xo + d/(2c))/y/A^. + (X + d/(2c)) 2 , and the prime is a shorthand for the deriva- 
tive w.r.t r. 

The first and simplest case we want to consider corresponds to the SDE (pi) with time inde- 
pendent parameter e > 0, D 2 > 0, and negative a. Whenever this conditions are satisfied, the 
process converges exponentially to the stationary regime with a typical relaxation time given 
by — l/o. Following [3] the stationary PDF can be computed in closed-form expression as 
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with v — 1 — 2a/ c, and the inverse tangent function continues smoothly at x > —2e/d. For 
illustrative purposes we fix the five free parameters as a = —20, b — d — e = 0.1, and c = 4.5, 
while the choice of g(t) and to is & t the moment irrelevant since we work directly with time 
t. As evident from equation (15) all moments of order n higher than or equal to v diverge: 
being v ~ 9.9, only the first nine lowest moments converge. In figure la we plot the time 
evolution of P(z t ,t) for increasing values of r = 0.01,0.05,0.1,1 as obtained by means of 

m = 2 13 , Az = —2z m i n /m, and At = 10 -3 ). For r = 1 we also plot 
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Figure 1: PDF of Z T for increasing values of r; solid line corresponds to the analytical stationary solution, 
dashed ones to FCA, while bars in Panel (a) and symbols in Panel (b) to MC simulation of the process for 
maturity r = 1. 



the histogram corresponding to MC simulation of the discrete process (parameter of the Euler 
scheme approximation: At = 10~ 3 , and number of MC paths Nuc = 10 6 ), while the solid 



line represents the analytical solution easily derived from equation (15). In figure lb we show 



the same results in log-linear scale to emphasize the tail region. The analytical information 
provides an overall check that the algorithm converges to the correct distribution, however far 
from stationary regime we have no precise information about the PDF shape. In [H] the scaling 
of the convergent moments is computed analytically, but it is known that the knowledge of the 
moments does not allow for a univocal reconstruction of the complete distribution. We can 
only rely on MC simulation, but sampling of low probability region requires on average a huge 
statistics (we need Nmc > 1/p to explore a p-probability region). At this point the advantages 
provided by the fast convolution based approach are evident, as clearly shown by both panels. 
FCA curve for r = 1 is in perfect agreement with the analytical prediction, both in central 
and tail regions. MC histograms agree as well, but these results are extremely noisy and very 
inaccurate for P(z T , r) < 10~ 4 . FCA based results are even more impressive looking at the 
computational time. Performances are strongly machine dependent, and for this reason we do 
not quote absolute times, but measured relative values: to obtain r = 1 bars MC takes ten 
times more than FCAPj As a consequence it needs 10 7 times more to reach the same accuracy 
at P(z T , r) ~ 1CT 10 level. Similar results are obtained for the slightly more complicated process 
used in [23J to model foreign exchange rate fluctuations. Their process is still mean reverting 
with a = -4.4 x 10 _1 , b = 0, c = 3.8 x 1(T 2 , and d = 3.04 x 10" 3 , though in this case the 
last parameter has a non trivially time dependence, e(r) = 6.08 x 10~ 5 + 6 x 10 _3 exp (-0.br). 
Lacking stationarity, every analytical information on the PDF is lost, yet we can see from 



figures 2a and 2b how the numerical PDF evolves with time and we verify a striking matching 
between MC and FCA results. Remarks similar to previous case apply. 

We now turn our attention to piecewise linear diffusion. The procedure is in this case a 
little bit subtle. While computation of the Lamperti transform of process (El) for H = 1/2 and 
integral time r = 2-y/t is still feasible providing qJ 
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3 Random numbers generators and FFT algorithms are provided by GNU Scientific Library. 
4 The sign function is defined according to the convention sign(0) = 0. 
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Figure 2: PDF of Z T for increasing values of r; lines correspond to FCA, while bars in Panel (a) and symbols 
in Panel (b) to MC simulations. 



the stochastic differential dZ T cannot be computed applying Ito Lemma straightforwardly. 
As a function of X T and r, Z T lacks necessary regularity condition for r = and X T = 0. 
However, both difficulties can be overcome. The X T process does not suffer any problem in 
r = 0, therefore we can evolve from X° to X 1 , and then exploit the one to one correspondence 
between X T and Z T . Moving from r = to r = Ar, X 1 remains delta distributed around zero, 
and the same holds true for Z 1 . For r > Ar the time derivative dZ T /dr needed in dZ T can 
be readily computed. The difficulty arising with the computation of d 2 Z T /dX 2 in zero can be 
dealt by replacing the absolute value with the smooth approximation 
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where in view of our application on a discrete grid the last approximation can be justified for 
sufficiently large k. We thus end with the following expression for the dynamics of Z T 



dZ T ~ sign(Z T ) 
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Numerical results concerning this last process are reported in figure 3a (linear scale) and in 

figure 

P{x,t) 



(log-linear scale). We study the dependence of P(z, r) on l/(cr 2 e). Indeed for |x| 3> 1, 
j exp [—2 |x| /(a 2 er)], and the value of the coefficient in the exponential function is 
crucial to asses the convergence of the expectation of exp(a;) with respect to P(x,r). We fix 
r = 1, a 2 — 1, and e = 0.5, 1,2. The leptokurtosis of the PDF increases as far as the value of 
e increases. Parameters for the Euler scheme approximation are fixed as in previous examples, 
while for the FCA we have slightly changed the value of Ar = 10~ 4 and m = 2 11 keeping 
^min = —10.24. For each one of the three cases we also plot the analytical prediction, since 
for piecewise linear processes the solution is known in closed form. Also in this last case the 
agreement between analytical and fast convolution PDF is totally satisfactory, while limitations 
of the MC approach are evident from symbols depicted in Panel (b). 
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Figure 3: PDF of Z T at time r = 1 for a 2 = 1 and e = 0.5,1,2. Panel (a): comparison between analytical 
expressions (solid lines) and FCA (dashed and dotted lines); Panel (b): comparison between MC histograms 
(symbols) and FCA. Log-linear curves have been shifted for readability. 

3 Financial applications 

In the second part of this work we present and discuss how results achieved in previous sessions 
can be exploited in finance, and in particular in the context of option pricing. For both quadratic 
and piecewise diffusion we briefly review how to set the correct risk neutral framework. Then, 
for explanatory purposes, we apply FCA to price European Plain Vanilla and geometric Asian 
options, but the approach can be extended to deal with different payoffs and different kind of 
boundary conditions. For the remaining of this paper, X t = In St — In St is the logarithmic 
return obtained from the stochastic process St describing the evolution of an asset price. As 
asset candidates we only consider equities and foreign exchange rates. 



3.1 Piecewise diffusion under risk-neutrality: Plain Vanilla pricing 

According to risk neutral valuation theory we need to find the dynamics of St or, equivalently, X t 
under the probability measure which makes all discounted asset prices martingales. Whenever 
the Novikov condition for the process under consideration is verified, Girsanov theorem gives 
the recipe for the equivalent measure, and it also explains how the dynamics of St coherently 
modifies. However, (32J [M] show how to compute the desired martingale directly from the 
Green function solving the FP equation associated to the dynamics 
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with S to = S . Just in the case of the original model of jU |3S], a delta hedged strategy allows 
to construct a locally risk neutral portfolio and to derive the partial differential equation 
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that is to be use to solve the pricing problem of a Plain Vanilla option O, with r risk free 
interest rate and for suitable boundary conditions. Introducing 0(S t ,t) = e r ( T ~^0(S t ,t) and 



substituting in equation (17), it is readily verified that the hat price satisfies an equation 



formally identical to the backward time FP equation associated with the dynamics (16) with 
ji — r. The fair price of a call option is therefore predicted to be 

/+oo 
dS T (S T - K) + G Q (S T ,T; S ,t ) = e^" 4 ^ [(S T - K) + \S ] , 
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where G^ is the Green function solving the FP equation in the risk neutral framework, and K 
is the strike price. The dynamics of X t under the new probability measure reads 
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At variance with equation (El), a non trivial drift term appears and some comments are manda- 
tory. As recognized by McCauley and collaborators, whenever the drift depends explicitly on 
X t there is no way to preserve scaling properties. However, in order to exploit the analytical 
information provided by equation ([5]), corresponding to the Green function G^(Xt, T; 0, 0) for 
the process Q, they provide arguable arguments in order to replace the drift with a constant. 
We are instead equipped with a computationally efficient algorithm, and so we can get rid of 
this approximation and price options directly with the process (18). As in sectional we switch 
to integral time 



dX T 



— -dr \X T dr + crW- + e\X 7 

2/2 2 ' ' V 2 ' 



Xn = 



and compute the Call option price as 

0(S Q ,to) = S°E Q [{e x ^ - e k )+|X ] = S°E Q [{e x ^m) 

sign(2 T(T) )i 



5, 



D 



T \ Z t{T) +V -2~ 



C 

r(T) 



k^+ 



Zo] 



m—l 



s °az y: 

j=0 



dz r(T ) 



— e 



p Q (2; r (T)|2;o) 



reAr 
2 



,nQ 



(19) 



with r(T) = 2a/T, discounted price Sq = e r( ^ T ^Sq an d log-moneyness k = ln(X/5o). The 
vector P™ 5 ^ of transition probability between zq and z" under the risk neutral measure Q has 



to be computed with the fast convolution procedure described in section 2.1 



3.2 Exotic options: the geometric Asian case 



Formula (19) can be extended to deal with payoffs with different dependence on S T , e.g. dig- 



ital options, covered call or strongly non linear function f(S T ), the only constraint being 
E® [/(5V) | So] < oo. The case of a functional payoff depending multiplicatively on the price 
along the path, i.e. /([S T ]) = Tl^of*^ 1 ), is just slightly more complicated but in fact can 
be easily managed, see [16] for an application to bond pricing. In this section we address the 
problem of pricing a geometric Asian option, which requires to compute the expected value 



E^ 



^T^/ t ln ^-ir) + |S 



(20) 



Being the positive part function non linear, previous expression is quite tricky to be evaluated 
and requires some careful manipulations. Defining r = 2(y/l — y/to), the Asian price is given 
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by 



OA(So,tc 



S^ 



gT-t JO V ■■! 



|-+V*oJX r /dr' _ k 



\Xn 



(21) 



We exploit the discretization of the t(T) time interval in n equally spaced intervals of amplitude 
At, and we replace the integral expression with a finite sum U n = Y^j=i (j'At/2 + v^o) & l n - 
We then introduce the ancillary variables {U 1 , . . . , U n } satisfying the following recursive relation 



U 



j+i 



i + l V 2 i + 1 



X 



i+l 



(22) 



for i = 1, . . . , n— 1 and f/ 1 = (4r + i/to) X 1 . Exploiting the one to one correspondence between 
X T and Z T , it is possible to rewrite the Asian price as 



O A (S ,t ) = S» / du n A(u n )p$(u n ) = S { 



du r 



dz n A(u n )p® z (u n ,z* 



with A(u n ) 



/T-VtQ , 



e k ) . The only unknown quantity in previous expression is the 



joint distribution of U n and Z n , whose computation requires a recursive relation allowing to 
propagate p^ z (u\ z l ) to p^ z (u t+1 , z t+1 ) with the associated initial time condition p u " z (u 1 , z 1 ) = 
5(z 1 )5(u 1 ). The following equation holds 



p^ z (u\z l+1 )= / dz l p$ z {u\z l )^{z l+l \z l ) 



(23) 



and to proceed it is useful to explicit the dependence of X l+1 on Z l+1 in equation (22) 

^ sign(Z i+1 )± 



U i+1 = ^U* 



2 ~ i+l 



\Z i + 1 \ 



(24) 



Z i+l = Z i+1 

where U l+l is coupled with the dummy variable Z l+1 . From previous relations we have 

d(u\ z l+1 ] 



pUu 1+ \z 1+1 ) 



d(u i+1 ,z i 



+i> 



^ Z («V + V +1 ),^ +1 ) 



(25) 



where the Jacobian is equal to (i+l)/i. Therefore starting from the distribution pjj Z (u 1 , z 1 ), and 
following the above procedure, after n — 1 steps we obtain the desired p^ z (u n , z n ). Introducing 
m^-node grid for Z % and m^-node grid for U\ we can approximate the distribution pjj z (u\ z l ) 
with a my x mz matrix P* fc , the row index j running over the nodes of U l , the column index 
k over those of Z % . The Asian price can therefore be approximated by 



mu—l 



mz-l 



O A (S ,t ) ~ S°AuAz J2 A K) J2 P 



)k 



3=0 



fc=0 



Since the time evolution corresponding to equation (23) is the most computationally intensive 
operation implicit in previous approximation, we can perform it at each node w* by means of 
FCA. The overall numerical complexity of the algorithm is essentially linear in the total number 
of grid nodes, i.e 0{n x mu x m z \og 2 m z ). 
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3.3 The Vellekoop-Nieuwenhuis-Borland model 



The geometric Asian case just described is useful also in view of the last application we present, 
which is related to the model for the stock price dynamics introduced in j6j[7]. Borland model 
tries to generalize the standard Black&Scholes to account for the empirical evidences of fat 
tailed return distributions, still keeping a closed form formula for the price of Plain Vanilla 
instruments. It is a sort of hybrid between a stochastic volatility model and the standard 
Black&Scholes: the volatility is stochastic, but the stock price and the volatility itself are 
driven by the same Brownian motion. Despite its theoretical elegance, the model contains 
some weaknesses extensively analyzed in [3H]. They raised two main objections by proving the 
Borland model to suffer arbitrage opportunities and diverging payoff expectation, as a conse- 
quence of the thickness of the tails. However, they preserved the main idea of the model and 
they amended it from both drawbacks. In their modified version they introduce the following 
dynamics 



dS t = fiStdt + aS t dVL t , S to = S 
dft t = E(ft t ,*)dW t , Q t0 = Q , 



(26) 
(27) 



with 



£(JM 



A-fP(ftt,t)' 





t> 
t = 0, 



and P(Q t ,t) = i F (l + aPtSlf) 



1 



(28) 



with ft = [(l-«)(2-«)t]"^, N t = A/y/K, A = yjT(±-§)/r(±), a e (0, §) and 
to > 0. In [38] the existence of a solution for equation (27) is proved, and it is shown how 
the unconditional distribution (i.e. Q to = for t = 0) reduces to the generalized Student-i 
distribution. In general the conditional distribution deviates from it. The log- return X t satisfies 
the equation 



X T = X t + fi(T -t)- ^a 2 J S 2 ( S , n s )ds + a(n T - tt t ) , 



for T > t > to- They verify that sufficient conditions hold for the applicability of the Girsanov 
theorem, and they derive the risk neutral dynamics 



dSt = rS t dt + aS t dnf , with d£l? = Z(n t ,t)dW® . 



(29) 



This last process does not suffer anymore of previous problems, however St does not satisfy the 
Markov property, yet only jointly with Q t . In addition, the price of Plain Vanilla instruments 
cannot be given in closed form formula. Indeed, according to pricing theory we have 



Oc(Sq, flo,t c 



S?E< 



r{T-t )+a(fi T -fl )- 



■ 2 f/tf(n 3 ,s)ds V 



Sq, fig 



and the expectation can only be computed via numerical techniques. 



Given (28), we observe that the equation governing the evolution of d£lf belongs to the 
class of quadratic diffusion processes Q through the identifications a = b = d = 0, c = 
a/ [(1 - a) (2 - a)], e(t) = [(1 - a) (2 - a)]^ t 2 ^ 2 ~ a \ and g(t) = t. Switching to the integral 
time t = hit /to for t Q > 0, and recalling equation (14), Z T is readily computed 



Z T 



[asinh (C ait0iT fi T ) - asinh (C a)t0iT fi )] , 



(30) 



12 



0.5 



FCA e = 0.5 
FCA e = 1.0 
FCA e = 2.0 




io- 



io- 





* MC e = 0.5 - 


^^^^ a ^^^ i ^. 


° MCe = 1.0 - 


O OO-. 


° MC e = 2.0 - 




\\ 



2.5 



(b) 



Figure 4: Piecewise diffusion: Risk neutral PDF of Z T at time r = 1 for r = 0.03, <r 2 = 1, and e = 0.5, 1,2. 
Comparison between Monte Carlo histograms (symbols) and FCA (dashed and dotted lines); in Panel (b) curves 
have been shifted for readability. 



where C a ^ hT = v/a/3t e T ^ 2 a \ The price of a Plain Vanilla instrument can be computed as 



O c (S Q ,Z ,to) = S£E 



) r(T-to)+ CT [n(Z T(T) )-n(Z )]-^f^[e(T)- e (to)]-£# / T(T) n(Z r ,) 2 dr' _ k 



<So, Zq 
(31) 



with t(T) = InT — lnt - Defining the set of ancillary variables {U 1 , . . . , U n } satisfing the 
recursive relation 



U l+1 = U l + AtVL(Z 



i+l\2 



(32) 



with U l = Arfi(Z 1 ) 2 , the formal analogy with the Asian case discussed in previous section is 



evident. Computation of the expectation in (31) requires estimation of the joint probability 



Pu Z (u n , z n ). The procedure is identical to the Asian case; equation (23) is still valid, while the 
system (24) has to be coherently modified in 



jji+i = U i + 



7i+l yi+1 



At 
r 2 

a,t Q ,(i + l)AT 



sinh 2 [y/cZ t+1 + asinh (C atto ^ i+1)AT fi(Z ))] 



The Jacobian in equation (25) simplifies to one, and, eventually, we can approximate the Plain 
Vanilla price as 

mjj — 1 mz — 1 

O c (S ,Z ,t )^S°AuAz JT ]T C( M ",4)P} fc Q , 

where C(u],z%) = h^-M^^ I^s a/2 ds + ain(z-)-n(z o) y 2(1 .^ ( 2 2 . Q) »" _ A ^ and compute it by 
means of fast convolution. 



3.4 Numerical results 

In this final section we sum up numerical results for the financial applications described in 



paragraphs 3.1, 3.2, and 3.3 
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IN) 

a. 




Figure 5: Piecewise diffusion: Bivariate risk neutral PDF of Z T and U T , and their corresponding marginals; 
e = 2, a 2 = 1, to = 0, and r = 1. In Panel (b) comparison between Fast Convolution PDF of Z and MC 
histogram. 



Whenever we switch to the risk neutral measure for the piecewise linear process, corrections 
terms in the SDE appear and an analytical expression for the density is not available anymore. 
Numerical simulation is mandatory, and the FCA algorithm, being both faster and much more 
efficient, is a natural competitor to MC approach. The transformed Z T process is enriched by 
the risk neutral correction (the last but one term in squared brackets) 



dZ T ~ sign(Z T ) 




where r is the risk free rate. In figures 4a and 4b we draw risk neutral PDFs for r = 0.03 and 



remaining parameters as in figures 3a and 3b The effect of the additional terms is evident from 
their comparison. In particular, it is remarkable the increase of the skewness induced by the 
risk neutral correction from linear plots in Panel (a). Turning our attention to the pricing of 

plots the joint density p(j z (u, z), and associated marginals for t = 0, 
2. Parameters of the fast convolution are z min = —10.24, uiz = 2 10 , 

3 



5a 



Asian options, Figure 
1, a 2 = 1, and e = 
= -2.56, mi, = 2 11 



r 



u 



and At = 10 . In figure 5b we compare the marginal PDF of Z T 



and analogously to paragraph 2J2 the agreement between FCA and MC is striking in the central 
region. 

As far as the pricing under Vellekoop-Nieuwenhuis-Borland model is concerned, we start 



plotting in figure 
2 io 



6a 



the joint bivariate density p^ z (u,z) for parameter values z min = —10.24, 
-5.12, m u = 2 11 , At = 10" 3 ', fin = 0, a = 0.1, U = 0.2, and T = 0.7. 



m-z = * , u min = — o.iz, mu — ^ , <-*i : - j-u , js — u, u. — u.i, l 
We notice that fast convolution algorithm correctly predicts a non negative support for the U T 
variable, even though the numerical grid spans uniformly the interval [u-mm , —u m \-n ]- In figure 6b 
we compare the distribution of Z T obtained by means of FCA and MC, finding perfect matching, 
and we also plot Q PDF, easily derived given the relationship between the two variables, see 
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Figure 6: Vellekoop-Nieuwenhuis-Borland model: Bivariate risk neutral PDF of Z and I/, and their correspond- 
ing marginals; a = 0.1, Oo = 0, to = 0.2, and T = 0.7. In Panel (b) plot of the fast convolution PDFs of Z and 
£1 and MC histogram of Z. 



equation (30). In light of the agreement between the two numerical procedures, we can use 



the FCA approach to efficiently price European Call options, as explained in paragraph 3.3 



In this respect in figures 7a| 7b 8a and 8b we present our results in terms of implied Black- 
Scholes volatilities. Our choices of the parameters are S to = 100, r = 0.03, a = 0.3, t = 0.2, 
tt = 0, 0.5, a = 0.1, 0.4, K e [70, 130], and T-t G [0.5, 2]. MC bands at 95% Confidence Level 
are plotted as dashed lines for the shortest time to maturity, T — 1 = 0.5 with Nmc = 5 x 10 7 . 
FCA and MC volatility curves are fully consistent. As expected surfaces exhibit a volatility 
smile, more pronounced for small maturities and for Q values deviating from zero. As already 
pointed out by Vellekoop and Niueuwenhuis, a wider variety of volatility surfaces and flexibility 
of the model can be obtained by playing with different values of Qq. 



4 Conclusions and perspectives 

In this paper we have addressed the problem of investigating performances of the fast convo- 
lution algorithm introduced by [20]. Choosing different specifications of the stochastic process, 
this has been carried out both with the reconstruction of conditional probability densities at 
different time horizons and with the computation of prices of financial derivatives. FCA is an 
efficient grid algorithm relying on restating functional integrals as sequences of ordinary finite 
dimensional integrals, and on converting the stochastic process to a unitary diffusion one by 
means of the Lamperti transform. A bright formulation of the problem, then, allows those 
integrals to be evaluated efficiently by the use of fast Fourier transform techniques. 

The stochastic processes we have investigated belongs to two classes of multiplicative noise 
processes: the family of quadratic diffusion, see [HI [19] , and piecewise linear diffusions, see [321 
[TJ . The analysis performed in this work provides a natural complement to the analytical results 
obtained in [TTJ, where closed form solutions for the stationary PDF and for the convergent 
moments at arbitrary time had been obtained. We have detailed a step by step numerical 
procedure able to provide an accurate estimate for the probability distribution of the process 
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Figure 7: FCA implied volatility surfaces, a = 0.1, t = 0.2, Panel (a) Qo = 0, and Panel (b) Oq = 0.5; dashed 
lines for T — io = 0.5 correspond to 95% Confidence Level from MC simulation. 




35 £n 



Figure 8: FCA implied volatility surfaces, a = 0.4, to = 0.2, Panel (a) Hq = 0, and Panel (b) 17o = 0.5; dashed 
lines for T — t = 0.5 correspond to 95% Confidence Level from MC simulation. 
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even far from the stationary regime. Similar results have been found for the piecewise diffusion. 
In this latter case, if the dynamics is enriched with a non trivial drift term, scaling properties 
are not preserved anymore and every analytical information is lost. Being this exactly the 
situation we faced when switching to the risk neutral setting, FCA proved to be a very efficient 
and reliable approach to the problem of option pricing. A detailed empirical analysis for 
different specifications of the parameter values documents the superiority of the FCA approach 
to standard Monte Carlo simulations. We have also demonstrated the flexibility of the approach 
when dealing with exotic instruments, and exploited the formal analogy between geometric 
Asian option pricing and Plain Vanilla pricing under Vellekoop-Nieuwenhuis-Borland dynamics. 
Being an interesting hybrid between a geometric Brownian motion and a stochastic volatility 
model, the latter provides a realistic description of the dynamics implied in the option market. 
FCA is able to numerically reproduce a rich variety of implied volatility surfaces improving the 
standard Monte Carlo approach. 

Since, as documented, FCA turns out to be highly successful also in the case of the Vellekoop- 
Nieuwenhuis-Borland model, a natural perspective is to concentrate future research efforts on 
the extension of FCA to higher dimensional stochastic systems. This is precisely the case of 
continuous time stochastic volatility models, see [21] . These models provide a flexible framework 
when modeling volatility, and they allow to reproduce several observed statistical regularities. 
For this reason they are nowadays extensively exploited by quantitative sectors of banks and 
financial institutions. Given the ability of the fast convolution to reconstruct densities over tail 
regions, and of the investigated models to generate leptokurtic and scaling distributions, the 
present approach is naturally suited for application in the context of financial risk management, 
e.g. Value-at-Risk and coherent risk measures computation, see [271 1361 ITUl 19]. 
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